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Abstract 

Compressed sensing typically deals with the estimation of a system input from its noise-corrupted linear mea- 
surements, where the number of measurements is smaller than the number of input components. The performance 
of the estimation process is usually quantified by some standard error metric such as squared error or support set 
error. In this correspondence, we consider a noisy compressed sensing problem with any arbitrary error metric. We 
propose a simple, fast, and highly general algorithm that estimates the original signal by minimizing the error metric 
defined by the user. We verify that our algorithm is optimal owing to the decoupling principle, and we describe a 
general method to compute the fundamental information-theoretic performance Umit for any error metric. We provide 
two example metrics — minimum mean absolute error and minimum mean support error — and give the theoretical 
performance limits for these two cases. Experimental results show that our algorithm outperforms methods such 
as relaxed belief propagation (relaxed BP) and compressive sampling matching pursuit (CoSaMP), and reaches the 
suggested theoretical limits for our two example metrics. 

Index Terms 

BeUef propagation, compressed sensing, error metric, estimation theory. 

I. Introduction 

A. Motivation 

Compressed sensing [1-3] has received a great deal of attention in recent years, because it deals with signal 
reconstruction problems with far fewer samples than required by the Nyquist rate, greatly reducing the sampUng 
rates required in signal processing appUcations, such as cameras, medical scanners, and high speed radar [4]. 

Consider a linear system, 

w = (1) 

where the system input X € is independent and identically distributed (i.i.d.), and the random linear mixing 
matrix $ G j^mxat known (typically M < N). The vector w £ is called the measurement of x, and is 
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passed through a bank of separable channels characterized by conditional distributions, 

M 

/Y|w(y|w) = Yl fY\w{yi\wi)- (2) 

i=l 

Note that the channels are general and are not restricted to Gaussian. We observe the channel output y, and want 
to estimate the original input signal x from y and 

The performance of the estimation process is often characterized by some error metric that quantifies the distance 
between the estimated and the original signals. For a signal x and its estimate x, both of length N, the error between 
them is the summation over the component-wise errors, 



For example, if the metric is absolute error, then d{xi, Xi) = \xi — Xi\; for squared error, d{xi, Xi) = {xi — Xi)^. 

Squared error is one of the most popular error metrics in various problems, due to many of its mathematical 
advantages. For example, minimum mean squared error (MMSE) estimation provides both variance and bias 
information about an estimator [5], and in the Gaussian case it is Unear and thus often easy to implement [6]. 
However, there are applications where MMSE estimation is inappropriate, for example because it is sensitive to 
outliers [7,8]. Therefore, alternative error metrics, such as mean absolute error (median), mean cubic error, or 
Hamming distance, are used instead. Considering the significance of various types of error metrics other than 
squared error, a general estimation algorithm that can minimize any desired error metric is of interest. 

B. Related work 

As mentioned above, squared error is most conmionly used as the error metric in estimation problems given 
by (1) and (2). Mean-square optimal analysis and algorithms were introduced in [9-13] to estimate a signal from 
measurements corrupted by Gaussian noise; in [14—16], further discussions were made about the circumstances 
where the output channel is arbitrary, while, again, the MMSE estimator was put forth. Another line of work, 
based on a greedy algorithm called orthogonal matching pursuit, was presented in [17, 18] where the mean squared 
error decreases over iterations. Absolute error is also under intense study in signal estimation. For example, an 
efficient sparse recovery scheme that minimizes the absolute error was provided in [19, 20]; in [21], a fundamental 
analysis was offered on the minimal number of measurements required while keeping the estimation error within 
a certain range, and absolute error was one of the metrics concerned. Support recovery error is another metric 
of great importance, for example because it relates to properties of the measurement matrices [22]. The authors 
of [22-24] discussed the support error rate when recovering a sparse signal from its noisy measurements; support- 
related performance metrics were appUed in the derivations of theoretical limits on the sampUng rate for signal 
recovery [25, 26]. The readers may notice that previous work only paid attention to limited types of error metrics. 
What if absolute error, cubic error, or other non-standard metrics are required in a certain application? 
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Fig. 1. Tanner graph for relaxed belief propagation. 
C. Contributions 

In this correspondence: (/) we suggest a Bayesian estimation algorithm that minimizes an arbitrary error metric; 
(//) we prove that the algorithm is optimal; (///) we study the fundamental information-theoretic performance limit 
of an estimation for a given metric; and (iv) we derive the performance limits for minimum mean absolute error and 
minimum mean support error estimators. This algorithm applies the relaxed belief propagation (BP) method [10, 
15] and the decoupling principle [9, 11, 12, 14]. It is simple and fast, and it reconstructs the original signal based on 
minimizing the expected error metric of interest. This is convenient for users who desire to recover the original signal 
with a non-standard error metric. Simulation results show that our algorithm outperforms algorithms such as relaxed 
BP [15], which is optimal for squared error, and compressive sampUng matching pursuit (CoSaMP) [18], a greedy 
reconstruction algorithm. Moreover, we compare our algorithm with the suggested theoretical limits for minimum 
mean absolute error (MMAE) and minimum mean support error (MMSuE), and illustrate that our algorithm is 
optimal. 

The remainder of the correspondence is arranged as follows: we review the idea of relaxed BP and the decoupling 
principle in Section II, and then describe our estimation algorithm and discuss its performance in Section III. 
Simulation results are given in Section IV, while conclusions appear in Section V. Some mathematical details 
appear in appendices. 

II. Review of Relaxed Belief Propagation 

Before describing the estimation algorithm, a review of the relaxed BP method [10, 15] is helpful. 

Belief Propagation (BP) [27] is an iterative method used to compute the marginals of a Bayesian network. 
Consider the bipartite graph, called a Tanner or factor graph, shown in Figure 1, where circles represent random 
variables (called variable nodes), and related variables are connected through functions (represented by squares, 
called/acfor nodes or function nodes) that indicate their dependence. In standard BP, there are two types of messages 
passed through the nodes: messages from variable nodes to factor nodes, rrix^y, and messages from factor nodes 
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to variable nodes, rUy^x- If we denote the set of function nodes connected to the variable x by N{x), the set of 
variable nodes connected to the function y by N{y), and the factor function at node y by ^y, then the two types 
of messages are defined as follows [27]: 

keN{x)\y 

n^v^x = 5Z ^yme^y. 

£eN{y)\x 

Inspired by the basic BP idea described above, the authors of [10, 14] developed iterative algorithms for estimation 
problems in linear mixing systems. In the Tanner graph, an input vector x = [xi,X2, ■■■,xn] is associated with 
the variable nodes (input nodes), and the output vector y = [yi, 2/21 •■•) 2/m] is associated with the function nodes 
(output nodes). If 7^ 0, then nodes xj and y, are cormected to an edge where the set of such edges E is 
defined as E = : ^ 0}. 

In standard BP methods [28,29], the distribution functions of Xj and ui-, as well as the channel distribution 
function fY\w{yi\'^i) were set to be the messages passed along the graph, but it is difficult to compute those 
distributions, making the standard BP method computationally expensive. In [10], a simplified algorithm, called 
relaxed belief propagation, was suggested. In this algorithm, means and variances replace the distribution functions 
themselves and serve as the messages passed through nodes in the Tanner graph, greatly reducing the computational 
complexity. In [14—16,30], this method was extended to a more general case where the channel is not necessarily 
Gaussian. 

An important result for BP in compressed sensing is that, as the signal dimension — > 00 and M/N is fixed, 
the BP process decouples the linear mixing channels (2) to a series of parallel scalar Gaussian channels [9, 11, 12, 
14]. In other words, estimating the input signal x from the output y of a linear mixing system (1), (2) is statistically 
equivalent to estimating each input entry Xj from a scalar Gaussian channel output qf 

qj=Xj+Vj, for i = 1, 2, 3, AT, (4) 

where each channel's additive Gaussian noise vj is jV(0,/u) distributed, and n satisfies Tanaka's fixed point 
equation [9,11,12,16,28,31]. 

III. Estimation Algorithm 

A. Algorithm 

Figure 2 illustrates the structure of our metric-optimal estimation algorithm (dashed box). The inputs of the 
algorithm are: (;) a distribution function /x(x), the prior of the original input x; (/;) a vector q — {qi,q2, ...,qN), 
the outputs of the scalar Gaussian channels computed by relaxed BP [15]; (///) a scalar /x, the variance of the 
Gaussian noise in (4); and (iv) an error metric function D{x,x) specified by the user. 
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Fig. 2. The structure of the metric-optimal estimation algorithm. 

Now that we know that the scalar channels have additive Gaussian noise, and that the variances of the noise are 
all /i, we can compute the conditional probability density function /x|Q(x|q) from Bayes' rule: 

„ , I . /Qix(q|x)/x(x) 

^^1^^"'^^ = /Q(q) 

/Qix(q|x)/x(x) 

//Q|x(q|x)/x(x)dx' 

Given an error metric D{% x), the optimal estimand Xopt is generated by minimizing the conditional expectation 
of the error metric E[D{x,x)\q\, which is easy to compute using /x|Q(x|q): 

£[£>(x,x)|q] = y£)(x,x)/x|Q(x|q)rfx 



/ 



£)(x,x) 



V(27rM) 



exp 



l|q-x|| 

2iJ, 



dx. 



Then, 



Xopt = argmJn J D{x,x)fy^iQ{x\q)dx. 



(6) 



(7) 



We solve for Xopt based on characteristics of the metric function D{x, x), and we differentiate between continuous 
and discrete cases. 

1) D{x,x) is continuous: Under a condition that D{x,x) is continuous and differentiable almost everywhere, 
we simply take the derivative of (6), set it equal to zero, and solve for Xop^, 

a/D(x,x)/x|Q(x|q)(ix 



0. 



(8) 



Because both the error metric function D{x, x) and the conditional probability /x|Q(x|q) are separable, the problem 
reduces to scalar estimation [6]. The derivative (8) is solved in a component- wise fashion: 

d J d{xi,Xi)f{xi\qi)dx. 



dxi 



= 0, 



(9) 



for each / G {1, 2, 3, N}. We discuss absolute error as an example of this case in Section 111-C. 

2) D{x, x) is discrete: If the error metric function is discrete, then we evaluate the minimizer Xopt by exhaustive 
search. We discuss support recovery error as an example of this case in Section III-D. 
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B. Theoretical results 

Having discussed the algorithm, we now give a theoretical justification for its performance. Note that the 
decoupling principle that we use in our proofs is based on the repUca method in statistical physics, and is not 
rigorous [9]. Therefore, our results are given as claims. On the other hand, repUca analyses in compressed sensing 
can be made rigorous in some cases; see for example Bayati and Montanari [32]. 

Claim 1: Given the system model described by (1), (2) and an error metric £'(x,x) of the form defined by (3), 
as the signal dimension N ^ oo and the measurement ratio M/N is fixed, the optimal estimand of the input signal 
is given by 

Xopt = arg niin E [D{% x) |q] , (10) 

X 

where the vector q is the output of the decoupled Gaussian scalar channel (4). 

Proof: From the main result by Guo and Verdu [9], linear mixed channels can be decoupled to a bank of scalar 
Gaussian channels qi — Xi + Vi where Vi ^ A/'(0,/i) for i E {1,2,...,A^} as iV oo. Then, ,fxi\Y{xi\y) gives 
the same distribution as fx-\Q-{xi\qi). In other words, once we know the value of /i, estimating each Xi from all 
channel outputs y = (yi,y2, vm) is equivalent to estimating Xi from the corresponding scalar channel output qi. 
Therefore, an estimator based on minimizing the conditional expectation of the error metric, -E (£'(x, x)|q), gives 
the best possible result. ■ 
Claim 1 states that, in the large system limit, the estimator satisfying (10) or (7) is optimal, because it minimizes 
the conditional expectation of the error metric. The key point in the estimation problem is to obtain the posterior 
fxi\'Y- Fortunately, the relaxed BP algorithm provides an asymptotically optimal method to decouple the mixing 
channels and thus an equivalent posterior fxi\Qi{xi\qi) can be computed easily, and our algorithm utilizes this 
convenient feature. 

Following Claim 1, we can compute the minimum expected error achievable by any estimation algorithm for any 
error metric _D(xopt,x). This minimum expected error is the fundamental information-theoretic performance limit 
of interest in this problem; no estimation algorithm can out-perform this limit. At the same time, we will see for 
two example error metrics that our BP-based algorithm matches the performance of the information-theoretic limit, 
and is thus optimal. 

Claim 2: With the optimal estimand Xopt determined by (7), the minimum mean user-defined error (MMUE) is 
given by 

MMUEl/x,.) = /^^^^ (7=S°P (-^)) /-W^ 

where R{-) represents the range of a variable, and ji is the variance of the noise of the scalar Gaussian channel 
(4). 
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Proof: 



MMUE(/x,m) =^^[^(xopt,x)] 



Ed 

R(Q) 



£;[£'(xopt,x)|q] 



/Q(q)c^q 



£;[i:»(xopt,x)|q]/Q(q)dq 

R{Q) 



£'(xopt,x)/x|Q(x|q)rfx /Q(q)rfq 

fl(Q) \JRm ) 

./Q|x(q|x)/x(x) \ f I 
, , D(Xopt,x) — ' rfx /Q(q)rfq 

_R(Q) \^fi(x) jQiqj y 



D(xopt,x)/Q|x(q|x)/x(x)dx dq 

i?(Q) V"'«(x) y 

■ 

Using both claims, we further analyze the estimation performance limits for two well-known error metrics in 
Sections IH-C and IH-D. 



C. Absolute error 

Because the MMSE is the mean of the conditional distribution, the outliers in the set of data may corrupt the 
estimation, and in this case the minimum mean absolute error (MMAE) is a good alternative. For absolute error, 
d^{xi,Xi) = |fi — a;i|, and we have the following corollary describing the performance limit of an MMAE 
estimator, where the proof is given in Appendix A. 

Corollary 1: For a system modeled by (1), (2), the minimum mean absolute error estimator achieves 

/ + 00 / /-Si, MMAE ^- + 00 \ 

(/ {-Xi)fxi\Qi{xi\li)dXi+ I Xifx,\Q,{xi\qi)dxi\ fQ,{qi)dqi, (11) 

where Xi (respectively, qi) is the input (respectively, output) of the scalar Gaussian channel (4), Ji,MMAE satisfies 
■/^,™ae ^Xi\Qi{xi\qi)dxi = \, and fxi\Qi{xi\qi) is a function of /x, following (5). 

D. Support Recovery Error 

In some applications in compressed sensing, correctly estimating the locations where the data has non-zero values 
is almost as important as estimating the exact values of the data; it is a standard model selection error criterion [22]. 
The process of estimating the non-zero locations is called support recovery. Support recovery error is defined as 
follows, and this metric function is discrete, 

c^support — Pr (xor(a;j, Xj)) , 
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where 

0, if a;,: = and Xi = 



xor(xi,Xi) 



0, if a;j 7^ and x » ^ 

1, if = and Xi^O 
1, if Xi 7^ and Xi — 

Corollary 2: For a system modeled by (1), (2), where /x is an i.i.d. sparse Gaussian prior such that Pr(Xi 
0) = p and ^ ~ AA(0, a^), the minimum mean support error estimator achieves 

MMSuE(/x, ,)=N.il-p). erfc (y^) + AT • p • erf {^^^) , (12) 

where 

crV/i y P y 

Corollary 2 is proved in Appendix B. 

IV. Numerical Results 

Some numerical results are shown in this section to illustrate the performance of our estimation algorithm 
when minimizing a user-defined error metric. The Matlab implementation of our algorithm can be found at http: 
//people.engr.ncsu.edu/dzbaron/software/arb_metric/. 

We test our estimation algorithm on two linear systems modeled by (1) and (2): (;) Gaussian input and Gaussian 
channel; (ii) WeibuU input and Poisson channel. In both cases, the input's length N is 10,000, and its sparsity rate is 
3%, meaning that the entries of the input vector are non-zero with probabihty 3%, and zero otherwise. The matrix 
# we use is Bernoulli(0.5) distributed, and is normahzed to have unit-norm rows. In the first case, the non-zero 
input entries are Af(0, 1) distributed, and the Gaussian noise is J\f(0,3 ■ 10~*) distributed, i.e., the signal to noise 
ratio (SNR) is 20dB. In the second case, the non-zero input entries are WeibuU distributed. 



f{xi;X,k) 



Xi < 



where A = 1 and k = 0.5. The Poisson channel is 

fY\w{y>i) = ^ " ^ , foralUe {1,2,...,M}, 

where the scahng factor of the input is a = 100. 

In order to illustrate that our estimation algorithm is suitable for reasonable error metrics, we considered absolute 
error and two other non-standard metrics: 

N 

Errorp = ^|%-Xj|P, 

where p = 0.5 or 1.5. 
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1000 2000 3000 4000 5000 6000 7000 
Number of measurements M 

Fig. 3. Comparison of the metric-optimal estimation algorithm, relaxed BP, and CoSaMP. (Sparse Gaussian input and Gaussian 
channel; sparsity rate= 3%; input length N=10,000; SNR=20dB.) 

We compare our algorithm with the relaxed BP [15] and CoSaMP [18] algorithms. In Figure 3 and Figure 4, lines 
marked with "metric-optimal" present the errors of our estimation algorithm, and lines marked with "Relaxed BP" 
(respectively, "CoSaMP") show the errors of the relaxed BP (respectively, CoSaMP) algorithm. Each point in the 
figure is an average of 100 experiments with the same parameters. Because the Poisson channel is not an additive 
noise channel and is not suitable for CoSaMP, the "MAE" and the "Errori.s" lines for "CoSaMP" in Figure 4 
appear beyond the scope of vertical axis. It can be seen that our metric-optimal algorithm outperforms the other 
two methods. 

To demonstrate the theoretical analysis of our algorithm in Sections III-C and III-D, we compare our MMAE 
estimation results with the theoretical limit (11) in Figure 5, where the integrations are computed numerically. In 
Figure 6, we compare our MMSuE estimator with the theoretical limit (12), where the value of /i is acquired 
numerically from the relaxed BP method [33] with 20 iterations. In both figures, each point on the "metric-optimal" 
line is generated by averaging 40 experiments with the same parameters. It is shown from both figures that the two 
lines are on top of each other. Therefore our estimation algorithm reaches both limits and is optimal. 

V. Conclusion 

In this correspondence, we introduced a simple and fast estimation algorithm that deals with arbitrary error metrics 
in noisy compressed sensing. We verified that the algorithm is optimal in the large system limit, and provided a 
general method to compute the minimum expected error achievable by any estimation algorithm for a user-defined 
error metric. We started with the scalar Gaussian channel model of the relaxed BP algorithm and extended it to a 
method that is applicable to any user-defined error metric. We discussed two error metric examples, MMAE and 
MMSuE, and gave the theoretical performance limits for them. We also illustrated numerically that our algorithm 
reaches these two theoretical limits, and outperforms the relaxed BP and the CoSaMP methods. 
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7000 



Fig. 4. Comparison of the metric-optimal estimation algorithm, relaxed BP, and CoSaMP. The "MAE" and the "Errori.^ " lines for 
"CoSaMP" appear beyond the scope of vertical axis. (Sparse WeibuU input and Poisson channel; sparsity rate— 3%; input length 
N= 10,000; input scaling factor a = 100.) 



Appendix A 
Proof of Corollary 1 



When 



equations (7) and (8) solve for the MMAE estimand, xmmae- We must satisfy (9), 

d E{\xi - Xi\ \ Qi) 



dxi 



0, 



(13) 



for each i = {1,2, • • • ,7V}. But 



E{\x^ - Xi\\qi 



poo 

I Pr(|fi - Xi\ > t\qi)dt 
Jo 



Pt(x, -Xi> t\qi)dt + / Pr{x^ - < ~t\q,)dt 
Jo 



Piixi < ti\q,)dti + / Pr(xi > t2\qi)dt2, 

) J Xi 

where changes of variables ti = Xi — t and t2 —Xi+t are applied in (14). Using (13) and (14), we need 

Pv{xi < XiMMAE\qi) - PriXi > Xi^MMAE\qi) = 0, 

and thus XiMMAE is given as the median of the conditional statistics fxi\Qi{xi\qi), 

fXi\QAx^h)dxi = / fx^\Q^ix^\q^)dxi = -. 



(14) 
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Fig. 5. Comparison of the minimum mean absolute error (MMAE) estimator and the theoretical limit (1 1). The two lines are on 
top of each other. (Sparse Gaussian input and Gaussian channel; sparsity rate— 3%; input length N=10,000; SNR=20dB.) 

Then, the conditional mean absolute error is, 

E[\Xi^MMAE - Xi\\qi] 

+ 00 

|Sz,MMAE - Xi\fx^\Q^{Xi\qi)dXi 

oo 

^i,MMAE /' + 00 

{XiMMAE - Xi)fx,\Q,{^i\Qi)dXi + / {Xi- Xi^mMAe) fx,\Q,{XiWi)dXi 



'OO 

^i,MMAE 



2^ i, MMAE 



+00 



{-Xi)fxi\Qiixi\qi)dx, + / xJx^\Q,ixi\qi)dxi. 



Therefore, the MMAE for location i, MMAEi(/x, , m), is 
MMAE,(/x,, m) = ^^[|^i,MMAE - 

+ 00 

E[\xiMMAE ~Xi\\ qt]fQ,{qi)dqi 



— OO 

+ 00 / pXiMMAE 



— OO \ J —OG 



{-Xi)fx^\Q^{xi\qi)dxi 



+ 00 



Xtfx,\Q.{^i\li)'i^i fqMijdqi 



We note in passing that the integrations can be evaluated numerically in an implementation. 

Because the input x is i.i.d., and the decoupled scalar channels have the same parameter /i, the values of MMAE^ 
for all i E {1, 2, A^} are the same, and the overall MMAE is 

MMAE(/x, fi) = N- MMAE,{fx, , ^i). 



Appendix B 
Proof of Corollary 2 

Similar to the idea of giving a limit on support recovery error rate [23], we derive the MMSuE limit for the 
case where the input is real-valued and the matrix $ is rectangular [M < N). In the scalar Gaussian channel (4), 
we factor the sparse Gaussian input Xi into Xi — Ui ■ Bi, where Ui ^ A/'(0, cr^) and Bi ^ Bernoulli(p), i.e., 
Pr(i3j = 1) = p =1 — Pr(i?i = 0). The support recovery problem is the task of finding the maximum a-posteriori 
(MAP) estimation of B,. 
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1000 2000 3000 4000 5000 6000 7000 
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Fig. 6. Comparison of the minimum mean support error (MMSuE) estimator and the theoretical limit (12). The two lines are on 
top of each other. (Sparse Gaussian input and Gaussian channel; sparsity rate= 3%; input length N=10,000; SNR=20dB.) 



For our estimation algorithm, the conditional expectation of support recovery error is, 

Pr(B, = l|q,) if fe, = and 6, = 1 



E [d{xi,Xi)\qi] = < 



Pt{B, = 0\q,) if bi 1 and 6, = 

if 6^ = and 5, = 

if 6i = 1 and bi = 1 
The estimand 6j,opt minimizes E [d{xi,Xi)\qi\, which implies 

if Pr(S, = <Pr(B^ =0|(7,) 

1 if Pr(S, = >Pr(S, =0|q,) 

It is easy to see that /Q.|s.fe|0) -7V(0,/i) and A/'(0, cr^ + Then, 

Pr(i?i ^l\q^) 



/Q,|B,fe|l)Pr(i3, = l) 



and similarly 



/Q,|B.fe|l)Pr(^. = l) 
E6.=o,i fQ,\BMh)Pr{B, = b,) 
1 

l + ^y'xVMTTexp(-f .|^) 
/Q^|B.fa|0)Pr(i3, ^0) 

/q,(*) 

i + i^y^v^exp(-f ■ 

Therefore, Pr(i3i = l\qi) > Pr(i?i = 0\qi) implies 



Pr(H, = 0|gO = 



qt >T = 2- 



and vice versa. We can rewrite (15) as. 




(15) 
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By averaging over the range of Qi, we get the overall MMSuE, 

MMSuE(/x,;,m) = N ■ E[d,^,^pp^n{x^,Xi)] 
= AT y E {di^support{xi,Xi)\qi) fQ.{qi)dqi 

= nI Vr{Bi = Q\qi)fQMi)dqi + N [ PriBi = l\qi)fQ^{qi)dqi 
= N ■ Pt{B, ^0,qf >t)+N ■ Pr(B, = 1, qf < t) 
= N ■ Pr(52 > ^ 0)PriBi = 0) + N ■ Pr{qf < t\B, = l)Pr{Bi = 1) 
= Ml-p)-erfc(y|:)+A.p.erf(y^^). 
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